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Abstract 

We consider the steady state equations for a compressible fluid. Since we wish to solve for a range of 
speeds we must consider the equations in conservation form. For transonic speeds these equations are of 
mixed type. Hence, the usual approach is to add time derivatives to the steady state equations and then 
march these equations in time. One then adds a time derivative of the density to the continuity equation, a 
derivative of the momentum to the momentum equation and a derivative of the total energy to the energy 
equation. This choice is dictated by the time consistent equations. However, since we are only interested 
in the steady state this is not necessary. Thus we shall consider the possibilty of adding a time derivative 
of the pressure to the continuity equation and similar modifications for the energy equation. This can then 
be generalized to adding combinations of time derivatives to each equation since these vanish in the steady 
state. When using acceleration techniques such as residual smoothing, multigrid, etc. these are applied to 
the pressure rather than the density. Hence, the code duplicates the behavior of the incompressible equations 
for low speeds. 


‘This research was supported by the National Aeronautics and Space Administration under NASA Contract No. NASl- 
19480 while the second author was in residence at the Institute for Computer Applications in Science and Engineering (ICASE), 
NASA Langley Research Center, Hampton, VA 23681-0001. 
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Introduction 


It is well known, that it is difficult to solve the compressible equations for low Mach numbers. For an 
explicit scheme this is easily seen by inspecting the time steps. For stability, the time step must be chosen 
inversely proportional to the largest eigenvalue of the system which, for slow flows, is approximately the 
speed of sound, c. However, other waves are convected at the fluid speed, u , which is much slower. Hence, 
these waves don’t change very much over a time step. Thousands of time steps may be required to reach a 
steady state. Should one try a multigrid acceleration one finds that the same disparity in wave speeds slows 
down the multigrid acceleration. With an implicit method an ADI factorization is generally used so that 
one can easily invert the implicit factors. The use of ADI introduces factorization errors which again slows 
down the convergence rate when there are wave speeds of very different magnitudes [7] . 

We consider systems of the form 

Wt 4 fx + 9y = 0. 

Our analysis will be based on the linearized equations so that the conservation form does not appear in the 
analysis though it does appear in the final numerical approximation. This system is now replaced by 

P -1 U>< + fx + 9y = o, 
or in linearized form 

P~ l w t 4 Aw x 4 Bw y = 0, (I) 

with A and B constant matrices. 

In order for this system to be equivalent to the original system, in the steady state, we demand that 
P* 1 have an inverse. This only need be true in the flow regime under consideration. We shall see later that 
frequently P is singular at stagnation points and also along sonic lines. Thus, the final preconditioner will 
be smoothed out in the vicinity of points where M=0 or M=l. 

Assuming the steady state has a unique solution, it does not matter which system we march to a steady 
state. We shall later see that for the finite difference approximations the steady state solutions are not 
necessarily the same and usually the preconditioned system leads to a better behaved steady state. 


The Incompressible Limit 


The time dependent two dimensional Euler equations can be written as 


Pt 4 up x -f vp y 4 pa 2 {u x 4 v y ) 

Px 

u t 4 uu x 4 vu y H 

P 

v t 4 uv x 4 vv y 4 — 
P 

S t 4 uS x 4 vS y 


0 

0 

0 

0 


( 2 ) 


The form of this system is unchanged if we nondimensionalize the equations. From now on we shall 
assume that are nondimensional quantities where the dimensional variables are nondimensionalized 

by Following [4] we define € = ^ = M* . If the fluid is isentropic then 


P - 




( 3 ) 


Hence, as e goes to zero the speed of sound, a, goes to infinity and so the first equation in (2) reduces to 
U X 4 Vy = 0. 
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It was pointed out in ([10], [11]) that these equations can be symmetrized by using & as the independent 
variable rather than dp . Hence, we define <j> by d<f) = For isentropic flow both p and a are functions only 

of the density and so using (3) this can be integrated explicitly. This gives <fi = . As the Mach number 

2 € 

goes to zero <j> tends to infinity and therefore, Gustafsson and Stoor [4] subtract a constant and define 
— 1 

* = 2=1. (4) 

2 e 

This amounts to specifying the constant in the integration of d<j> from dp. They then prove, using energy 
methods, that 

&Pin compressible 
a<t>x — 

ax 

Hence, <j> and all its derivatives behave as O(M) as M — ► 0. Since p —+ 1 and using the definition of d<f> this 
is equivalent to 

dp compressible * ^Pi ncompressible (5) 

We now consider how to construct a matrix artificial viscosity that will enable us to reach the incom- 
pressible limit. Consider 

P Wt + f x -b Qy “ h [(Ql’Wx)x d" (Q2^y)y] (6) 

We wish to find the dependence of P and Qi on the Mach number as Af — ► 0 so that we get the proper 
convergence. We therefore consider the isentropic equations based on w = v) see (4). This has the 

symmetric form 


a n 

<*12 

^13 

\ 

a 12 

<*22 

<*23 

Wx 

a n 

<*23 

<*33 / 

I 


C &ll &12 3 \ 

b\2 b 22 ^23 I Wy = h [(Q\W X ) X + (Q2U>y)y] 

b\l &23 633 f 

As M — *■ 0 ,ai 2 and b 13 = 0(1/ M) while d</> — O(M) while all other quantities are bounded. Hence, the 
leading terms in the first equation are all 0 ( 1 / M) while they are 0 ( 1 ) for the second and third equations. 
Multiplying the first equation by M and taking the limit we get u x + v y for the space derivatives on the left 
hand side. Using d<j> = 0(M),du = 0(1), dv = 0(1) we see that a necessary condition for convergence as 
M —+ 0 is that P” 1 , Qi have the form 


p-\Qi,Q 2 ' 


^(m) 

0(A) 0(1) 0(1) 

0(i) 0 ( 1 ) 0 ( 1 ) 


The artificial viscosity matrices Qi are related to the preconditioner. Consider the one dimensional 
equation 

u t + Pf x ~ (\Q\u x ) x 

Let A = |£. Since we are updating Pf we should have Q = PA. However, this is not in conservation form 
at the steady state. Instead we consider artificial viscosities of the form 

P- x u t + f x ~(P~'(\PA\u x ) x 


U t + Pf x ~ P{P-\\PA\U X ) X 
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This would be equivalent to the original form if P were constant. Instead we have terms like PjP ^i+i/2 
that appear. 

We note that the conditions on the matrix (P” 1 |PD|) C omp are not satisfied by the non-preconditioned 
Roe matrices. Furthermore, even reasonable preconditioners need not satisfy these conditions. Consider, for 
example, the one dimensional system 

P ~ l w t -F Aw x = h(Qw x ) x 

A reasonable choice is P" 1 = \A\ i.e. P = \A~ 1 \. In this case all the wave speeds of PA are ±1. Now 
Q = P-'\PA\ = \A\\\A\-'A\ = \A\~ 


( 0 (£) 0 ( 1 ) 0 ( 1 ) 
0 ( 1 ) 0(jf) 0 ( 1 ) 

V 0(1) 0(1) 0(1) 


i.e. Q is the nonpreconditioned Roe matrix which does not have the desired property. We therefore conclude 
that for an upwind difference scheme the Riemann solver should be based on the preconditioned system and 
not the original scheme. In [3] plots are shown to illustrate the greatly improved accuracy for low Mach 
number flows when the Riemann solver is based on the preconditioning. Characteristics in the boundary 
conditions these should be based on the characteristics of the modified system and not the physical system. 
Preconditioning is even more important when using multigrid than with an explicit scheme. With the original 
system the disparity of the eigenvalues greatly affects the smoothing rates of the slow components and so 
slows down the multigrid method, [6]. 

We conclude from the above remarks that the steady state solution of the preconditioned system may 
be different from that of the physical system. Thus, on the finite difference level the preconditioning can 
improve the accuracy as well as the convergence rate. 

Algorithm 


In terms of the primitive variables the preconditioning we consider is: 


/ 

a 2 

0 

0 

0 

\ 

aau 

~w 

1 

0 

0 


aav 

P 2 

0 

1 

0 


V 

0 

0 

0 

1 

/ 


( U 

a 

0 

0 

\ 


a 

u 

0 

0 


+ 

0 

0 

u 

0 



V o 

0 

0 

u 

) 


( V 

0 

a 

0 

\ 


0 

V 

0 

0 


+ 

a 

0 

V 

0 



l o 

0 

0 

V 

/ 


/ * \ 
du 
dv 

V dS 

I ) 

du 

dv 

dS I, 

& \ 

r pa 

du 

dv 

\dS 
x ' y 


The nonpreconditioned case corresponds to /? 2 = a 2 , a = 0. Let q = uuji -j- uu> 2 > then the eigenvalues of PD 
are given by 

do = q (double) 

d ± = \{(l-a+f/a>)± (8) 


(1 — a + /J 2 /a 2 ) 2 + 4(ui 2 + 



3 



For general curvilinear coordinates, in the ”i” direction <+\ = uj 2 = — a;,,. The time step is bounded by 
VOL 

<4 ' 

Our ultimate goal is to have a compressible code that solves the incompressible equations when the input 
Mach number is zero. So we wish to use variables that give us the same result as an incompressible code on 
all levels of the algorithm, e.g. flux computation, boundary conditions, acceleration techniques, etc. Hence, 
we choose as our basic variables 


( p' \ 


P' \ 

pu 

, Q = 

pu 

pv 


pv 

E> yl 


V n' j 




pu 

\ 



pv 

\ 

F = 


pu 2 + p' 


3 G — 


puv 




puv 




pv + p 




pH'u 

J 



pH'v 

) 


dWp_ [aF dG 
di ^ dx ^ dy 


where 

, , _ 

P — P Poo '*oo — « 

7- 1 

E' = c pP (T - Too) — {p - Poo) + ' 2 " J 

= P + Poo - hoop 

pH* — E f + p ~ E + p - h^p 


We subtract the constants to keep the quantities in scale. Density is now calculated from the pressure and 
total energy. Because the modified energy E l also contains the density we get a quadratic equation for 
the density. Choosing the positive square root guarantees that the density is always positive. The residual 
smoothing and multigrid are applied to p f and E ' rather than p and E. Thus, we duplicate the treatment of 
the variables in a pseudo-compressible incompressible code. 


P p = 1+ A 


/ 1-s 

—uB 2 
— vB 2 

V -b 4 


where h 
Si = 
#2 = 
b 4 = 


CpT : 


1 

w 

Si 


G+Zioo 

t FB* 

G+hoo 

uvB? 

G+hoo 
— uB 1 

G + hoo 

uvB 2 

G+hoo 

x +B? 

G+h^ 
— vBi 

G+hoo 

uBa 

G+hoo 

vBa 

G + hoo 
— B 4 

G+hoo 

G+hoo 

G+hoo 

_ u 2 +v 2 
— 2 

,A = 

(G + hoo)P' 
h 


■, In the appendix we derive this form of P 


p- 


1 

(7 - l)/» 




0 2 


B\H‘ + 


a(u 2 + v 2 ) 

J 2 
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We choose 0 2 = min {majc [ft(u 2 + t; 2 ), , a 2 } where 0 min should have the units of speed. The 
choice of 0m% n is discussed in the result section. In all cases a = min j^l, ♦ We can evaluate this 

efficiently by defining 5 — A * (dp - ud {p u )+ v ^{(*>)- dE j — $J\)h • Then 
dPnevu 

d(/?'u) rieit; — d(/?u) or j^ B 2 uS 
d(pv^ new — B 2 vS 

d^new — dEorig — B4S 

These equations are given for the nondimensionalized variables. The nondimensionalization affects the 
convergence. In some codes, p and p are fixed in the far field. This implies that the speed of sound, a, is 
also bounded. As the Mach number goes to zero the pressure remains of order 1 while the velocities go to 
zero. Alternatively, one can nondimensionalize so that the velocities are of order 1 in the far field and then 
the pressure and speed of sound go to infinity, unless one subtracts an appropriate constant, 

The boundary conditions at the far field boundary, for subsonic flow, are based on the one dimensional 
theory of characteristics in the direction normal to the boundary. The preconditioning changes the form of 
these characteristic variables. In differential form they are given by 

1 f 0 2 

Ri = du ~Wv {l ~ Q ~^ ) 

j dp 

R2 = du — 2^2 ( u(1 — a — “2 ) 

+^(«( l-a+^)) 2 + 4 (l-^)/? 2 ^ dp 1 

where u is the component of the velocity normal to the boundary. If we consider low Mach numbers then 
we can approximate these by 

dp' dp' 

R l= du + —, R 2 = du 

P0 P0 

which is the same as for the incompressible case. Hence, at inflow R\, v (tangential velocity) and S are 
specified while R 2 is extrapolated from the interior. We then calculate u (normal velocity) and the pressure 
from R\ and R 2 and then the density and total energy. At outflow the role of specified and extrapolated 
quantities is reversed. At solid boundaries the normal momentum equation is used which is not affected by 
the preconditioning. 

Computational Results 

The solution is advanced by a explicit Runge-Kutta method ([5], [8]) with residual smoothing and multigrid 
and no enthalpy damping. In all cases three levels of FMG multigrid were used with 50 Runge-Kutta cycles 
on the coarser grids. Hence, all plots show the convergence for two sets of 50 cycles and then the convergence 
on the finest mesh. The plots are of the convergence rate of the residual of the continuity equation. For the 
original code this was updated for the density while in the preconditioned code it is updated for the pressure. 
Nevertheless, in the steady state the residual of the continuity equation should be the same except for the 
change in the artificial viscosity between the two algorithms. All cases were run with a matrix viscosity. 

We first present two calculations for inviscid flow about a NACA 0012. We use a 224 x 32 C mesh 
and three levels of multigrid. The first calculation is for inflow conditions M = 0.01, a = 1.25° . In this 
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case we see that the residual asymptotes without the use of preconditioning and that the preconditioning 
dramatically increases the rate of convergence. The use of the preconditioning adds only a few percent 
to the total computational time. In the second case we consider the same geometry but with an inflow 
of M — 0.7, a = 1.25°. We have also done M = 0.8, a = 1.25° which results in a minor slowing of the 
convergence rate. The preconditioned residual is the dotted line and the original code is the solid line. 
Different parameters for the time step and residual smoothing are needed with and without preconditioning. 
For inviscid cases we can choose /? m * n as zero while for the viscous cases /? min = 0 Ay/u^ + For the 
transonic cases the lift and drag coefficients are changed only minimally by the preconditioning. 

We next consider viscous flow about a RAE2822 airfoil on a 320 x 64 C mesh and 5 levels of multigrid on 
the finest grid with = 0.01, a = 2.79° using a Baldwin-Lomax turbulence model with Re — 6.5 million. 
The residual history is presented in figure 3. Again the standard code converges very slowly for these low 
Mach numbers. In figure 4 we present both the preconditioned residual (dashed line) and the original code 
(solid line) for the same case but Moo = 0.73. For viscous cases we choose 0 m in = 0.4. Again, for the 
transonic cases the lift and drag are changed by about 2 percent by the preconditioning. For the very low 
Mach numbers the lift and drag coefficients never converged for the non-preconditioned algorithm and seem 
to have significant errors. The preconditioned code gives much better agrees for lift and drag for low Mach 
numbers. 

We conclude with a three dimensional case, inviscid flow about an ONERA wing. In figure 5 we display 
the convergence rate for the continuity equation (normalized by the initial residual) for Mach numbers .10, 
.05 and .01. We see that the convergence rate is independent of the inflow Mach number. In figure 6 we 
plot the lift coefficient for the same case. We again see that the lift coefficient is essentially independent 
of the Mach number except for some slight compressibility effects. Without preconditioning there are large 
variations in the lift for this set of Mach numbers. 
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Appendix 


To find P p we begin with the preconditioner P 5 for the variables dWs = du, dv, dS)*, with dS = 
dp — a 2 dp. We then transform to dW s = (dp, dti, dv , d 5 ) t by multiplying all elements in the first row of the 
matrix by pa and every element in the first column by This gives 


p-1 _ 


/ 

1 j* 

aau 

P 2 

aav 

P 2 

V 0 


0 0 0 \ 

1 0 0 
0 1 0 
001/ 


Pc = 


( 


V 


ii 

a 2 

— av, 
a 

— av 
a 

0 


0 0 0 \ 

1 0 0 
0 1 0 
001/ 


We then transform to the conservation variables W c = (p, pu, pv, E ) x . This is given by dW c — T\dWs • Let 

r = 7- 1 


/ - 1 


T\ = 


0 
0 

P 

\ ^ pu pv 


U 


0 

p 

0 


\ 


Ml I 

2 / 


/ TG 


TC l = 


\ TG 


-Tu -Tv r \ 
A 0 0 

0 ' / 

-r« -rr r / 


0 


where G = - and a is the speed of sound. This gives the preconditioner in conservation variables. 


^ _ ( 7 - 1 ) 0 0 2 -o!) 


_ ( 7 -l)(^ : 

2 -(l+or)a 2 


a* 

) W 2 


a 4 

conservation variables P c = 

TiP 

st; 1 . 



/ l + GQi 

-uQ 1 

-vQi 

Qx \ 

p 

«gq 2 

1 — u 2 

£?2 

-uvQ 2 

uQ 2 


| vGQ 2 

-*a^Q2 

1 - t; 2 Q 2 

vQ 2 


V gr 

— uR 

—vR 

1 + -ft / 

We next change 

& 

0 

3 

6 

c = (p,pu,p, £) variables 

dW' c = T 2 dW c . 






f 1 0 

0 1 \ 





0 1 

0 0 




t 2 = 

0 0 

1 0 





^ ^00 0 

0 1 ) 





( 1 + GQi 

-uQi 

-vQi 

-Qi \ 

p' — 

UGQ2 

1 — u 2 

o 2 

-uvQ 2 

uQ 2 

* c ~ 

vGQ 2 

— 

1 - v 7 Q 2 

vQ 2 


\ RG 

—uR 

—vR 

l + R / 


Then P' c = T 2 (T 1 P S TJ 1 )T2 1 
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We finally change only the time derivatives to the variables W p by dW p — T 3 dW' c . 


T a = 


( T{G + h x ) 
0 
0 

V 0 


-r« 

l 

o 

o 


-r\> 

0 

1 

o 


r \ 
o 

0 

1 / 


7T 1 = 


/ 


n G+&„) 

o 

o 

o 


G+h, 

1 

0 

0 


G +h 0 
0 
1 
0 


G+fco 

0 


To summarize, we begin with 

„-i dW s , „ dW s , 0 dW s „ 
Pc +As— fBs— =0 


dt dx 

and transform to 


dy 


Pr .w + w ! + fl w L = 0 








We then transform , in conservation form, to the prime variables where E is replaced by E f = E — /iooP- 


Finally we then have that P p 1 = P' 
Thus, 


- 1 rp— 1 


r 2 TiPs~ 1 T 1 - 1 T 2 - 1 r 3 - 1 . or P p = T 3 P' c _1 = T 3 T 2 Ti P sTf 1 T 2 - 1 


P p = 


/ A 
-uAf?2 
— i / A £?2 
V -A5 4 


1 + ^ 2 b 2 


h 

u £b. 


h 

UV0 2 B 7 
h 2 

1 | tflgga 

^ 2 * 
V0 2 B< 

h 


-U0 2 B2 

K 

-v0 2 B 7 


where all quantities were defined in the text. 

We next show how to convert any preconditioner given in streamline coordinates and (&, du, dv, dS) 
coordinates to conservative variables in Cartesian (not streamwise) coordinates. We shall do this in two 
dimensions but the extension to three dimensions is straightforward. Assume we are given a preconditioner 
in streamline coordinates and du, dv ) dS) coordinates Ps given by 


Pq = 


U 


u~ l = 



< Pn 

Pl2 

0 

0 

\ 


P 21 

p 22 

0 

0 



0 

0 

P 33 

0 



< 0 

0 

0 

Pss j 

define rotation matrices U, 

( i 

0 

0 

0 \ 



0 cosO 

sinO 

0 



o - 

sinO 

cosO 

0 



0 

0 

0 

1 } 



1 

0 

0 

0 \ 


0 

cosd 

— sinO 0 


0 

sind 

cosO 

0 


0 

0 

0 

1 

) 


Let q 2 — ir + v 2 . To get the streamwise direction we shall choose 


cosO — 


yju 2 + V 2 


sinO = 


vV 4- v 2 
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Then the preconditioner in Cartesian coordinates is given by Pear = U *PsU and 

/ Pn Pia“ Pu~ q 0 \ 

Pi if (Pii-P 33 m 0 

P 2 l^ ( P 2 2 ^33 ) P 2 2^+P33p 0 

0 0 0 P 55 / 

^ <2l] Q 12 Ql3 0 ^ 

_ Qn Q 22 Ql 3 0 

Q 31 Q 32 Q 33 0 

^ 0 0 0 Q 55 J 

We next introduce conservative variables W c as given in the appendix by the transformation T\. The 
preconditioner for conservative variables is then given by P c = TiPcarT'i 1 
We now define the following quantities 


Vi = 

y 2 - 

y 3 = 

Y a = 

n = 
g = 

L = 
Z\\ = 

Z 12 = 

2\z = 

Z 14 — 

^15 = 

^51 = 


Z 52 = 

^53 = 

Z54 = 

^55 = 


tiQ 2 i T i^Q 3 i + u;Q4i 
w Q 2 2 + ^Q 32 + ^$42 
uQ 23 + ^Q 33 + W<243 
uQ 24 + ^^34 + w £?44 

(7 - i)(Qs 5 - Qn) 

(7-l)g 2 

2 

(7 “ l)p 


r ^2 


1 

pa 2 [ 
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55 


1 

uY 5 

+ 

Q\2 
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P 

1 

■vY s 


Q 13 

— 0 
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pa 2 



p 

1 

0 

tnYs 

+ 

Q 14 


pa- p 
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J/Z n +GYi -^2-^3 


u?y4 


+<?55( J 


7 - 1 


) 


h Z 12 - (7 - i)p«yi + y 2 - UQ 55 
HZ13 - (7 - l)pwVl + ^3 - uQss 
H — (7 — l)ptuYi + V 4 — nQ 55 
HZ 15 + (7 - IW 55 
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Then 


Pc 


/ ^11 %12 

uZn -|- GQ21 — ^2 U Z \2 ~ RuQ 21 4 * Q22 

vZu + <jQ33 ~ *3 v %12 ~ RUQ 3 1 + Q32 

w Zll -h GQ4 1 ~ V4 t^Zi2 — ^ u Q 41 4- Q42 

\ ^51 ^52 


^13 

uZ!3 — HVQ22 + Q 23 
vZj 3 - HVQ 31 + Q 33 
t^Z!3 - RvQu + Q 43 
^53 


^14 

UZ14 - Htx;Q21 + Q24 
vZ i4 - RwQ 31 -f Q34 
wZii ~ R^Qa + Q44 
^54 


^15 

xiZis + RQ 21 
vZ]5 + .RQ 31 
4- RQ41 
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Figure 3: Convergence rate for viscous flow about a RAE2822 airfoil with M 0 0 = 0.01 and a = 2.79° Solid 
line is original algorithm and dashed line is the preconditioned scheme 



Figure 4: Convergence rate for viscous flow about a RAE2822 airfoil with M oo = 0.73 and a = 2.79°. Solid 
line is original algorithm and dashed line is the preconditioned scheme 
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Figure 5: Convergence rate for inviscid flow about ONERA wing, M 0 0 = .10, .05, .01, a = 3.06° 
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Figure 6: Lift coefficient for inviscid flow about a ONERA wing. 
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